

<!DOCTYPE html>
<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
<head>
  <meta charset="utf-8">
  
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  
  <title>Sparse Linear Equations &mdash; CVXOPT User&#39;s Guide</title>
  

  
  
  
  

  
  <script type="text/javascript" src="_static/js/modernizr.min.js"></script>
  
    
      <script type="text/javascript" id="documentation_options" data-url_root="./" src="_static/documentation_options.js"></script>
        <script src="_static/jquery.js"></script>
        <script src="_static/underscore.js"></script>
        <script src="_static/doctools.js"></script>
        <script src="_static/language_data.js"></script>
    
    <script type="text/javascript" src="_static/js/theme.js"></script>

    

  
  <link rel="stylesheet" href="_static/css/theme.css" type="text/css" />
  <link rel="stylesheet" href="_static/pygments.css" type="text/css" />
    <link rel="search" title="Search" href="search.html" />
    <link rel="copyright" title="Copyright" href="copyright.html" />
    <link rel="next" title="Cone Programming" href="coneprog.html" />
    <link rel="prev" title="Discrete Transforms" href="fftw.html" /> 
</head>

<body class="wy-body-for-nav">

   
  <div class="wy-grid-for-nav">
    
    <nav data-toggle="wy-nav-shift" class="wy-nav-side">
      <div class="wy-side-scroll">
        <div class="wy-side-nav-search" >
          

          
            <a href="index.html" class="icon icon-home"> CVXOPT User's Guide
          

          
          </a>

          
            
            
              <div class="version">
                1.2.5
              </div>
            
          

          
<div role="search">
  <form id="rtd-search-form" class="wy-form" action="search.html" method="get">
    <input type="text" name="q" placeholder="Search docs" />
    <input type="hidden" name="check_keywords" value="yes" />
    <input type="hidden" name="area" value="default" />
  </form>
</div>

          
        </div>

        <div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">
          
            
            
              
            
            
              <ul class="current">
<li class="toctree-l1"><a class="reference internal" href="copyright.html">Copyright and License</a></li>
<li class="toctree-l1"><a class="reference internal" href="intro.html">Introduction</a></li>
<li class="toctree-l1"><a class="reference internal" href="matrices.html">Dense and Sparse Matrices</a></li>
<li class="toctree-l1"><a class="reference internal" href="blas.html">The BLAS Interface</a></li>
<li class="toctree-l1"><a class="reference internal" href="lapack.html">The LAPACK Interface</a></li>
<li class="toctree-l1"><a class="reference internal" href="fftw.html">Discrete Transforms</a></li>
<li class="toctree-l1 current"><a class="current reference internal" href="#">Sparse Linear Equations</a><ul>
<li class="toctree-l2"><a class="reference internal" href="#matrix-orderings">Matrix Orderings</a></li>
<li class="toctree-l2"><a class="reference internal" href="#general-linear-equations">General Linear Equations</a></li>
<li class="toctree-l2"><a class="reference internal" href="#positive-definite-linear-equations">Positive Definite Linear Equations</a></li>
<li class="toctree-l2"><a class="reference internal" href="#example-covariance-selection">Example: Covariance Selection</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="coneprog.html">Cone Programming</a></li>
<li class="toctree-l1"><a class="reference internal" href="solvers.html">Nonlinear Convex Optimization</a></li>
<li class="toctree-l1"><a class="reference internal" href="modeling.html">Modeling</a></li>
<li class="toctree-l1"><a class="reference internal" href="c-api.html">C API</a></li>
<li class="toctree-l1"><a class="reference internal" href="printing.html">Matrix Formatting</a></li>
<li class="toctree-l1"><a class="reference external" href="http://cvxopt.org">cvxopt.org</a></li>
</ul>

            
          
        </div>
      </div>
    </nav>

    <section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">

      
      <nav class="wy-nav-top" aria-label="top navigation">
        
          <i data-toggle="wy-nav-top" class="fa fa-bars"></i>
          <a href="index.html">CVXOPT User's Guide</a>
        
      </nav>


      <div class="wy-nav-content">
        
        <div class="rst-content">
        
          















<div role="navigation" aria-label="breadcrumbs navigation">

  <ul class="wy-breadcrumbs">
    
      <li><a href="index.html">Docs</a> &raquo;</li>
        
      <li>Sparse Linear Equations</li>
    
    
      <li class="wy-breadcrumbs-aside">
        
            
        
      </li>
    
  </ul>

  
  <hr/>
</div>
          <div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
           <div itemprop="articleBody">
            
  <div class="section" id="sparse-linear-equations">
<span id="c-spsolvers"></span><h1>Sparse Linear Equations<a class="headerlink" href="#sparse-linear-equations" title="Permalink to this headline">¶</a></h1>
<p>In this section we describe routines for solving sparse sets of linear
equations.</p>
<p>A real symmetric or complex Hermitian sparse matrix is stored as an
<a class="reference internal" href="matrices.html#cvxopt.spmatrix" title="cvxopt.spmatrix"><code class="xref py py-class docutils literal notranslate"><span class="pre">spmatrix</span></code></a> object <code class="docutils literal notranslate"><span class="pre">X</span></code>  of size
(<img class="math" src="_images/math/5a939c5280da7202ca4531f175a7780ad5e1f80a.png" alt="n"/>, <img class="math" src="_images/math/5a939c5280da7202ca4531f175a7780ad5e1f80a.png" alt="n"/>) and an
additional character argument <code class="docutils literal notranslate"><span class="pre">uplo</span></code> with possible values <code class="xref py py-const docutils literal notranslate"><span class="pre">'L'</span></code>
and <code class="xref py py-const docutils literal notranslate"><span class="pre">'U'</span></code>.  If <code class="docutils literal notranslate"><span class="pre">uplo</span></code> is <code class="xref py py-const docutils literal notranslate"><span class="pre">'L'</span></code>, the lower triangular part
of <code class="docutils literal notranslate"><span class="pre">X</span></code> contains the lower triangular part of the symmetric or Hermitian
matrix, and the upper triangular matrix of <code class="docutils literal notranslate"><span class="pre">X</span></code> is ignored.  If <code class="docutils literal notranslate"><span class="pre">uplo</span></code>
is <code class="xref py py-const docutils literal notranslate"><span class="pre">'U'</span></code>, the upper triangular part of <code class="docutils literal notranslate"><span class="pre">X</span></code> contains the upper
triangular part of the matrix, and the lower triangular matrix of <code class="docutils literal notranslate"><span class="pre">X</span></code> is
ignored.</p>
<p>A general sparse square matrix of order <img class="math" src="_images/math/5a939c5280da7202ca4531f175a7780ad5e1f80a.png" alt="n"/> is represented by an
<code class="xref py py-class docutils literal notranslate"><span class="pre">spmatrix</span></code> object of size (<img class="math" src="_images/math/5a939c5280da7202ca4531f175a7780ad5e1f80a.png" alt="n"/>, <img class="math" src="_images/math/5a939c5280da7202ca4531f175a7780ad5e1f80a.png" alt="n"/>).</p>
<p>Dense matrices, which appear as right-hand sides of equations, are
stored using the same conventions as in the BLAS and LAPACK modules.</p>
<div class="section" id="matrix-orderings">
<span id="s-orderings"></span><h2>Matrix Orderings<a class="headerlink" href="#matrix-orderings" title="Permalink to this headline">¶</a></h2>
<p>CVXOPT includes an interface to the AMD library for computing approximate
minimum degree orderings of sparse matrices.</p>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<ul class="simple">
<li><p>P. R. Amestoy, T. A. Davis, I. S. Duff,  Algorithm 837: AMD, An
Approximate Minimum Degree Ordering Algorithm, ACM Transactions on
Mathematical Software, 30(3), 381-388, 2004.</p></li>
</ul>
</div>
<dl class="py function">
<dt id="cvxopt.amd.order">
<code class="sig-prename descclassname">cvxopt.amd.</code><code class="sig-name descname">order</code><span class="sig-paren">(</span><em class="sig-param">A</em><span class="optional">[</span>, <em class="sig-param">uplo = 'L'</em><span class="optional">]</span><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.amd.order" title="Permalink to this definition">¶</a></dt>
<dd><p>Computes the approximate mimimum degree ordering of a symmetric  sparse
matrix <img class="math" src="_images/math/211284f68205c3e66773eaf026f32a0acdd3dfb3.png" alt="A"/>.  The ordering is returned as an integer dense matrix
with length equal to the order of <img class="math" src="_images/math/211284f68205c3e66773eaf026f32a0acdd3dfb3.png" alt="A"/>.  Its entries specify a
permutation that reduces fill-in during the Cholesky factorization.
More precisely, if <code class="docutils literal notranslate"><span class="pre">p</span> <span class="pre">=</span> <span class="pre">order(A)</span></code> , then <code class="docutils literal notranslate"><span class="pre">A[p,</span> <span class="pre">p]</span></code> has
sparser Cholesky factors than <code class="docutils literal notranslate"><span class="pre">A</span></code>.</p>
</dd></dl>

<p>As an example we consider the matrix</p>
<div class="math">
<p><img src="_images/math/d66b4ab25c0eef60d649ef489104789f43f45a0b.png" alt="\left[ \begin{array}{rrrr}
 10 &amp;  0 &amp; 3 &amp;  0 \\
  0 &amp;  5 &amp; 0 &amp; -2 \\
  3 &amp;  0 &amp; 5 &amp;  0 \\
  0 &amp; -2 &amp; 0 &amp;  2
\end{array}\right]."/></p>
</div><div class="doctest highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">spmatrix</span><span class="p">,</span> <span class="n">amd</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">A</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">([</span><span class="mi">10</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="p">,</span><span class="o">-</span><span class="mi">2</span><span class="p">,</span><span class="mi">5</span><span class="p">,</span><span class="mi">2</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">P</span> <span class="o">=</span> <span class="n">amd</span><span class="o">.</span><span class="n">order</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="nb">print</span><span class="p">(</span><span class="n">P</span><span class="p">)</span>
<span class="go">[ 1]</span>
<span class="go">[ 0]</span>
<span class="go">[ 2]</span>
<span class="go">[ 3]</span>
</pre></div>
</div>
</div>
<div class="section" id="general-linear-equations">
<span id="s-umfpack"></span><h2>General Linear Equations<a class="headerlink" href="#general-linear-equations" title="Permalink to this headline">¶</a></h2>
<p>The module <code class="xref py py-mod docutils literal notranslate"><span class="pre">cvxopt.umfpack</span></code> includes four functions for solving
sparse non-symmetric sets of linear equations.  They call routines from
the UMFPACK library, with all control options set to the default values
described in the UMFPACK user guide.</p>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<ul class="simple">
<li><p>T. A. Davis, Algorithm 832: UMFPACK – an unsymmetric-pattern
multifrontal method with a column pre-ordering strategy, ACM
Transactions on Mathematical Software, 30(2), 196-199, 2004.</p></li>
</ul>
</div>
<dl class="py function">
<dt id="cvxopt.umfpack.linsolve">
<code class="sig-prename descclassname">cvxopt.umfpack.</code><code class="sig-name descname">linsolve</code><span class="sig-paren">(</span><em class="sig-param">A</em>, <em class="sig-param">B</em><span class="optional">[</span>, <em class="sig-param">trans = 'N'</em><span class="optional">]</span><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.umfpack.linsolve" title="Permalink to this definition">¶</a></dt>
<dd><p>Solves a sparse set of linear equations</p>
<div class="math">
<p><img src="_images/math/908a083a0806ab219190adac99ec60da04c65f48.png" alt="AX &amp; = B \quad (\mathrm{trans} = \mathrm{'N'}), \\
A^TX &amp; = B \quad (\mathrm{trans} = \mathrm{'T'}), \\
A^HX &amp; = B \quad (\mathrm{trans} = \mathrm{'C'}),"/></p>
</div><p>where <img class="math" src="_images/math/211284f68205c3e66773eaf026f32a0acdd3dfb3.png" alt="A"/> is a sparse matrix and <img class="math" src="_images/math/4bc3e94a67870b41b7c20179693e889251e2c136.png" alt="B"/> is a dense matrix.
The arguments <code class="docutils literal notranslate"><span class="pre">A</span></code> and <code class="docutils literal notranslate"><span class="pre">B</span></code> must have the same type
(<code class="xref py py-const docutils literal notranslate"><span class="pre">'d'</span></code> or <code class="xref py py-const docutils literal notranslate"><span class="pre">'z'</span></code>) as <code class="docutils literal notranslate"><span class="pre">A</span></code>.  On exit <code class="docutils literal notranslate"><span class="pre">B</span></code> contains
the solution.  Raises an <code class="xref py py-exc docutils literal notranslate"><span class="pre">ArithmeticError</span></code> if the coefficient
matrix is singular.</p>
</dd></dl>

<p>In the following example we solve an equation with coefficient matrix</p>
<div class="math" id="equation-e-sp-adef">
<p><span class="eqno">(1)<a class="headerlink" href="#equation-e-sp-adef" title="Permalink to this equation">¶</a></span><img src="_images/math/5046b770f82749051a01f8f102c8e5ff69c4a3e4.png" alt="A = \left[\begin{array}{rrrrr}
    2 &amp; 3 &amp; 0 &amp; 0 &amp; 0 \\
    3 &amp; 0 &amp; 4 &amp; 0 &amp; 6 \\
    0 &amp;-1 &amp;-3 &amp; 2 &amp; 0 \\
    0 &amp; 0 &amp; 1 &amp; 0 &amp; 0 \\
    0 &amp; 4 &amp; 2 &amp; 0 &amp; 1
    \end{array}\right]."/></p>
</div><div class="doctest highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">spmatrix</span><span class="p">,</span> <span class="n">matrix</span><span class="p">,</span> <span class="n">umfpack</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">V</span> <span class="o">=</span> <span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">6</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">I</span> <span class="o">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span>  <span class="mi">2</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span>  <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">4</span><span class="p">]</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">J</span> <span class="o">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span>  <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span>  <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">]</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">A</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">(</span><span class="n">V</span><span class="p">,</span><span class="n">I</span><span class="p">,</span><span class="n">J</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">B</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">1.0</span><span class="p">,</span> <span class="p">(</span><span class="mi">5</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">umfpack</span><span class="o">.</span><span class="n">linsolve</span><span class="p">(</span><span class="n">A</span><span class="p">,</span><span class="n">B</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="nb">print</span><span class="p">(</span><span class="n">B</span><span class="p">)</span>
<span class="go">[ 5.79e-01]</span>
<span class="go">[-5.26e-02]</span>
<span class="go">[ 1.00e+00]</span>
<span class="go">[ 1.97e+00]</span>
<span class="go">[-7.89e-01]</span>
</pre></div>
</div>
<p>The function <a class="reference internal" href="#cvxopt.umfpack.linsolve" title="cvxopt.umfpack.linsolve"><code class="xref py py-func docutils literal notranslate"><span class="pre">linsolve</span></code></a>  is
equivalent to the following three functions called in sequence.</p>
<dl class="py function">
<dt id="cvxopt.umfpack.symbolic">
<code class="sig-prename descclassname">cvxopt.umfpack.</code><code class="sig-name descname">symbolic</code><span class="sig-paren">(</span><em class="sig-param"><span class="n">A</span></em><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.umfpack.symbolic" title="Permalink to this definition">¶</a></dt>
<dd><p>Reorders the columns of <code class="docutils literal notranslate"><span class="pre">A</span></code> to reduce fill-in and performs a symbolic
LU factorization.  <code class="docutils literal notranslate"><span class="pre">A</span></code> is a sparse, possibly rectangular, matrix.
Returns the symbolic factorization as an opaque C object that can be
passed on to <a class="reference internal" href="#cvxopt.umfpack.numeric" title="cvxopt.umfpack.numeric"><code class="xref py py-func docutils literal notranslate"><span class="pre">numeric</span></code></a>.</p>
</dd></dl>

<dl class="py function">
<dt id="cvxopt.umfpack.numeric">
<code class="sig-prename descclassname">cvxopt.umfpack.</code><code class="sig-name descname">numeric</code><span class="sig-paren">(</span><em class="sig-param"><span class="n">A</span></em>, <em class="sig-param"><span class="n">F</span></em><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.umfpack.numeric" title="Permalink to this definition">¶</a></dt>
<dd><p>Performs a numeric LU factorization of a sparse, possibly rectangular,
matrix <code class="docutils literal notranslate"><span class="pre">A</span></code>.   The argument <code class="docutils literal notranslate"><span class="pre">F</span></code> is the symbolic factorization
computed by <a class="reference internal" href="#cvxopt.umfpack.symbolic" title="cvxopt.umfpack.symbolic"><code class="xref py py-func docutils literal notranslate"><span class="pre">symbolic</span></code></a>
applied to the matrix <code class="docutils literal notranslate"><span class="pre">A</span></code>,
or another sparse matrix with the same sparsity pattern, dimensions,
and type.  The numeric factorization is returned as an opaque C object
that that can be passed on to
<a class="reference internal" href="#cvxopt.umfpack.solve" title="cvxopt.umfpack.solve"><code class="xref py py-func docutils literal notranslate"><span class="pre">solve</span></code></a>.  Raises an
<code class="xref py py-exc docutils literal notranslate"><span class="pre">ArithmeticError</span></code> if the matrix is singular.</p>
</dd></dl>

<dl class="py function">
<dt id="cvxopt.umfpack.solve">
<code class="sig-prename descclassname">cvxopt.umfpack.</code><code class="sig-name descname">solve</code><span class="sig-paren">(</span><em class="sig-param">A</em>, <em class="sig-param">F</em>, <em class="sig-param">B</em><span class="optional">[</span>, <em class="sig-param">trans = 'N'</em><span class="optional">]</span><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.umfpack.solve" title="Permalink to this definition">¶</a></dt>
<dd><p>Solves a set of linear equations</p>
<div class="math">
<p><img src="_images/math/908a083a0806ab219190adac99ec60da04c65f48.png" alt="AX &amp; = B \quad (\mathrm{trans} = \mathrm{'N'}), \\
A^TX &amp; = B \quad (\mathrm{trans} = \mathrm{'T'}), \\
A^HX &amp; = B \quad (\mathrm{trans} = \mathrm{'C'}),"/></p>
</div><p>where <img class="math" src="_images/math/211284f68205c3e66773eaf026f32a0acdd3dfb3.png" alt="A"/> is a sparse matrix and <img class="math" src="_images/math/4bc3e94a67870b41b7c20179693e889251e2c136.png" alt="B"/> is a dense matrix.
The arguments <code class="docutils literal notranslate"><span class="pre">A</span></code> and <code class="docutils literal notranslate"><span class="pre">B</span></code> must have the same type.  The argument
<code class="docutils literal notranslate"><span class="pre">F</span></code> is a numeric factorization computed
by <a class="reference internal" href="#cvxopt.umfpack.numeric" title="cvxopt.umfpack.numeric"><code class="xref py py-func docutils literal notranslate"><span class="pre">numeric</span></code></a>.
On exit <code class="docutils literal notranslate"><span class="pre">B</span></code> is overwritten by the
solution.</p>
</dd></dl>

<p>These separate functions are useful for solving several sets of linear
equations with the same coefficient matrix and different right-hand sides,
or with coefficient matrices that share the same sparsity pattern.
The symbolic factorization depends only on the sparsity pattern of
the matrix, and not on the numerical values of the nonzero coefficients.
The numerical factorization on the other hand depends on the sparsity
pattern of the matrix and on its the numerical values.</p>
<p>As an example, suppose <img class="math" src="_images/math/211284f68205c3e66773eaf026f32a0acdd3dfb3.png" alt="A"/> is the matrix <a class="reference internal" href="#equation-e-sp-adef">(1)</a> and</p>
<div class="math">
<p><img src="_images/math/7b694e55a0813fe78e01121684b12a40fad86d86.png" alt="B = \left[\begin{array}{rrrrr}
    4 &amp; 3 &amp; 0 &amp; 0 &amp; 0 \\
    3 &amp; 0 &amp; 4 &amp; 0 &amp; 6 \\
    0 &amp;-1 &amp;-3 &amp; 2 &amp; 0 \\
    0 &amp; 0 &amp; 1 &amp; 0 &amp; 0 \\
    0 &amp; 4 &amp; 2 &amp; 0 &amp; 2
    \end{array}\right],"/></p>
</div><p>which differs from <img class="math" src="_images/math/211284f68205c3e66773eaf026f32a0acdd3dfb3.png" alt="A"/> in its first and last entries.  The following
code computes</p>
<div class="math">
<p><img src="_images/math/2e7d885be2f8d9a7556b4e6dded7b42315803547.png" alt="\newcommand{\ones}{\mathbf 1}
x = A^{-T}B^{-1}A^{-1}\ones."/></p>
</div><div class="doctest highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">spmatrix</span><span class="p">,</span> <span class="n">matrix</span><span class="p">,</span> <span class="n">umfpack</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">VA</span> <span class="o">=</span> <span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">6</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">VB</span> <span class="o">=</span> <span class="p">[</span><span class="mi">4</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">6</span><span class="p">,</span> <span class="mi">2</span><span class="p">]</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">I</span> <span class="o">=</span>  <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span>  <span class="mi">2</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span>  <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">4</span><span class="p">]</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">J</span> <span class="o">=</span>  <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span>  <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span>  <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">]</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">A</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">(</span><span class="n">VA</span><span class="p">,</span> <span class="n">I</span><span class="p">,</span> <span class="n">J</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">B</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">(</span><span class="n">VB</span><span class="p">,</span> <span class="n">I</span><span class="p">,</span> <span class="n">J</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">x</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">1.0</span><span class="p">,</span> <span class="p">(</span><span class="mi">5</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">Fs</span> <span class="o">=</span> <span class="n">umfpack</span><span class="o">.</span><span class="n">symbolic</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">FA</span> <span class="o">=</span> <span class="n">umfpack</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">Fs</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">FB</span> <span class="o">=</span> <span class="n">umfpack</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">B</span><span class="p">,</span> <span class="n">Fs</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">umfpack</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">FA</span><span class="p">,</span> <span class="n">x</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">umfpack</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">B</span><span class="p">,</span> <span class="n">FB</span><span class="p">,</span> <span class="n">x</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">umfpack</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">FA</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">trans</span><span class="o">=</span><span class="s1">&#39;T&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="nb">print</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
<span class="go">[ 5.81e-01]</span>
<span class="go">[-2.37e-01]</span>
<span class="go">[ 1.63e+00]</span>
<span class="go">[ 8.07e+00]</span>
<span class="go">[-1.31e-01]</span>
</pre></div>
</div>
</div>
<div class="section" id="positive-definite-linear-equations">
<span id="s-cholmod"></span><h2>Positive Definite Linear Equations<a class="headerlink" href="#positive-definite-linear-equations" title="Permalink to this headline">¶</a></h2>
<p><code class="xref py py-mod docutils literal notranslate"><span class="pre">cvxopt.cholmod</span></code> is an interface to the Cholesky factorization routines
of the CHOLMOD package.  It includes functions for Cholesky factorization
of sparse positive definite matrices, and for solving sparse sets of linear
equations with positive definite matrices.
The routines can also be used for computing
<span class="raw-html">LDL<sup><small>T</small></sup></span>
(or
<span class="raw-html">LDL<sup><small>H</small></sup></span>
factorizations
of symmetric indefinite matrices (with <img class="math" src="_images/math/19eef1966f7c545af3ac8c0fa486974d873e3c65.png" alt="L"/> unit lower-triangular and
<img class="math" src="_images/math/0fcab9067b50b87e868c4fd70f213a086addb964.png" alt="D"/> diagonal and nonsingular) if such a factorization exists.</p>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<ul class="simple">
<li><p>Y. Chen, T. A. Davis, W. W. Hager, S. Rajamanickam,
Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization
and Update/Downdate, ACM Transactions on Mathematical Software,
35(3), 22:1-22:14, 2008.</p></li>
</ul>
</div>
<dl class="py function">
<dt id="cvxopt.cholmod.linsolve">
<code class="sig-prename descclassname">cvxopt.cholmod.</code><code class="sig-name descname">linsolve</code><span class="sig-paren">(</span><em class="sig-param">A</em>, <em class="sig-param">B</em><span class="optional">[</span>, <em class="sig-param">p = None</em>, <em class="sig-param">uplo = 'L'</em><span class="optional">]</span><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.cholmod.linsolve" title="Permalink to this definition">¶</a></dt>
<dd><p>Solves</p>
<div class="math">
<p><img src="_images/math/e44d29f0cbcc5fa7a0f3f83d3bddc988f9b022ea.png" alt="AX = B"/></p>
</div><p>with <img class="math" src="_images/math/211284f68205c3e66773eaf026f32a0acdd3dfb3.png" alt="A"/> sparse and real symmetric or complex Hermitian.</p>
<p><code class="docutils literal notranslate"><span class="pre">B</span></code> is a dense matrix of the same type as <code class="docutils literal notranslate"><span class="pre">A</span></code>.  On exit it
is overwritten with the solution.  The argument <code class="docutils literal notranslate"><span class="pre">p</span></code> is an integer
matrix with length equal to the order of <img class="math" src="_images/math/211284f68205c3e66773eaf026f32a0acdd3dfb3.png" alt="A"/>, and specifies an
optional reordering.
See the comment on
<code class="xref py py-attr docutils literal notranslate"><span class="pre">options['nmethods']</span></code> for details on which ordering is used
by CHOLMOD.</p>
<p>Raises an <code class="xref py py-exc docutils literal notranslate"><span class="pre">ArithmeticError</span></code> if the factorization does not exist.</p>
</dd></dl>

<p>As an  example, we solve</p>
<div class="math" id="equation-e-a-pd">
<p><span class="eqno">(2)<a class="headerlink" href="#equation-e-a-pd" title="Permalink to this equation">¶</a></span><img src="_images/math/aaf2039a0ca644d308dc578303aa096fd3775907.png" alt="\left[ \begin{array}{rrrr}
        10 &amp;  0 &amp; 3 &amp;  0 \\
         0 &amp;  5 &amp; 0 &amp; -2 \\
         3 &amp;  0 &amp; 5 &amp;  0 \\
         0 &amp; -2 &amp; 0 &amp;  2
    \end{array}\right] X =
    \left[ \begin{array}{cc}
         0 &amp; 4 \\ 1 &amp; 5 \\ 2 &amp; 6 \\ 3 &amp; 7
    \end{array} \right]."/></p>
</div><div class="doctest highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">matrix</span><span class="p">,</span> <span class="n">spmatrix</span><span class="p">,</span> <span class="n">cholmod</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">A</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">([</span><span class="mi">10</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">5</span><span class="p">,</span> <span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="mi">5</span><span class="p">,</span> <span class="mi">2</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">X</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="mi">8</span><span class="p">),</span> <span class="p">(</span><span class="mi">4</span><span class="p">,</span><span class="mi">2</span><span class="p">),</span> <span class="s1">&#39;d&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">cholmod</span><span class="o">.</span><span class="n">linsolve</span><span class="p">(</span><span class="n">A</span><span class="p">,</span><span class="n">X</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="nb">print</span><span class="p">(</span><span class="n">X</span><span class="p">)</span>
<span class="go">[-1.46e-01  4.88e-02]</span>
<span class="go">[ 1.33e+00  4.00e+00]</span>
<span class="go">[ 4.88e-01  1.17e+00]</span>
<span class="go">[ 2.83e+00  7.50e+00]</span>
</pre></div>
</div>
<dl class="py function">
<dt id="cvxopt.cholmod.splinsolve">
<code class="sig-prename descclassname">cvxopt.cholmod.</code><code class="sig-name descname">splinsolve</code><span class="sig-paren">(</span><em class="sig-param">A</em>, <em class="sig-param">B</em><span class="optional">[</span>, <em class="sig-param">p = None</em>, <em class="sig-param">uplo = 'L'</em><span class="optional">]</span><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.cholmod.splinsolve" title="Permalink to this definition">¶</a></dt>
<dd><p>Similar to <a class="reference internal" href="#cvxopt.cholmod.linsolve" title="cvxopt.cholmod.linsolve"><code class="xref py py-func docutils literal notranslate"><span class="pre">linsolve</span></code></a> except that
<code class="docutils literal notranslate"><span class="pre">B</span></code> is an <a class="reference internal" href="matrices.html#cvxopt.spmatrix" title="cvxopt.spmatrix"><code class="xref py py-class docutils literal notranslate"><span class="pre">spmatrix</span></code></a> and
that the solution is returned as an output argument (as a new
<code class="xref py py-class docutils literal notranslate"><span class="pre">spmatrix</span></code>).  <code class="docutils literal notranslate"><span class="pre">B</span></code> is not modified.
See the comment on
<code class="xref py py-attr docutils literal notranslate"><span class="pre">options['nmethods']</span></code> for details on which ordering is used
by CHOLMOD.</p>
</dd></dl>

<p>The following code computes the inverse of the coefficient matrix
in <a class="reference internal" href="#equation-e-a-pd">(2)</a> as a sparse matrix.</p>
<div class="doctest highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="n">X</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">splinsolve</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">spmatrix</span><span class="p">(</span><span class="mf">1.0</span><span class="p">,</span><span class="nb">range</span><span class="p">(</span><span class="mi">4</span><span class="p">),</span><span class="nb">range</span><span class="p">(</span><span class="mi">4</span><span class="p">)))</span>
<span class="gp">&gt;&gt;&gt; </span><span class="nb">print</span><span class="p">(</span><span class="n">X</span><span class="p">)</span>
<span class="go">[ 1.22e-01     0     -7.32e-02     0    ]</span>
<span class="go">[    0      3.33e-01     0      3.33e-01]</span>
<span class="go">[-7.32e-02     0      2.44e-01     0    ]</span>
<span class="go">[    0      3.33e-01     0      8.33e-01]</span>
</pre></div>
</div>
<p>The functions <a class="reference internal" href="#cvxopt.cholmod.linsolve" title="cvxopt.cholmod.linsolve"><code class="xref py py-func docutils literal notranslate"><span class="pre">linsolve</span></code></a> and
<a class="reference internal" href="#cvxopt.cholmod.splinsolve" title="cvxopt.cholmod.splinsolve"><code class="xref py py-func docutils literal notranslate"><span class="pre">splinsolve</span></code></a> are equivalent to
<a class="reference internal" href="#cvxopt.cholmod.symbolic" title="cvxopt.cholmod.symbolic"><code class="xref py py-func docutils literal notranslate"><span class="pre">symbolic</span></code></a> and
<a class="reference internal" href="#cvxopt.cholmod.numeric" title="cvxopt.cholmod.numeric"><code class="xref py py-func docutils literal notranslate"><span class="pre">numeric</span></code></a> called in sequence, followed by
<a class="reference internal" href="#cvxopt.cholmod.solve" title="cvxopt.cholmod.solve"><code class="xref py py-func docutils literal notranslate"><span class="pre">solve</span></code></a>, respectively,
<a class="reference internal" href="#cvxopt.cholmod.spsolve" title="cvxopt.cholmod.spsolve"><code class="xref py py-func docutils literal notranslate"><span class="pre">spsolve</span></code></a>.</p>
<dl class="py function">
<dt id="cvxopt.cholmod.symbolic">
<code class="sig-prename descclassname">cvxopt.cholmod.</code><code class="sig-name descname">symbolic</code><span class="sig-paren">(</span><em class="sig-param">A</em><span class="optional">[</span>, <em class="sig-param">p = None</em>, <em class="sig-param">uplo = 'L'</em><span class="optional">]</span><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.cholmod.symbolic" title="Permalink to this definition">¶</a></dt>
<dd><p>Performs a symbolic analysis of a sparse real symmetric or
complex Hermitian matrix <img class="math" src="_images/math/211284f68205c3e66773eaf026f32a0acdd3dfb3.png" alt="A"/> for one of the two factorizations:</p>
<div class="math" id="equation-e-chol-ll">
<p><span class="eqno">(3)<a class="headerlink" href="#equation-e-chol-ll" title="Permalink to this equation">¶</a></span><img src="_images/math/c9a3af023076f58fe1ee0b35f88cad484ba194fa.png" alt="PAP^T = LL^T, \qquad PAP^T = LL^H,"/></p>
</div><p>and</p>
<div class="math" id="equation-e-chol-ldl">
<p><span class="eqno">(4)<a class="headerlink" href="#equation-e-chol-ldl" title="Permalink to this equation">¶</a></span><img src="_images/math/340c54d7b4f301eb582e5ce306364e1cc4fe20b1.png" alt="PAP^T = LDL^T, \qquad PAP^T = LDL^H,"/></p>
</div><p>where <img class="math" src="_images/math/c2aa3dff9bffb099e9dff196fd36aed56ec16baf.png" alt="P"/> is a permutation matrix, <img class="math" src="_images/math/19eef1966f7c545af3ac8c0fa486974d873e3c65.png" alt="L"/> is lower triangular
(unit lower triangular in the second factorization), and <img class="math" src="_images/math/0fcab9067b50b87e868c4fd70f213a086addb964.png" alt="D"/> is
nonsingular diagonal.  The type of factorization depends on the value
of <code class="xref py py-attr docutils literal notranslate"><span class="pre">options['supernodal']</span></code> (see below).</p>
<p>If <code class="docutils literal notranslate"><span class="pre">uplo</span></code> is <code class="xref py py-const docutils literal notranslate"><span class="pre">'L'</span></code>, only the lower triangular part of <code class="docutils literal notranslate"><span class="pre">A</span></code>
is accessed and the upper triangular part is ignored.
If <code class="docutils literal notranslate"><span class="pre">uplo</span></code> is <code class="xref py py-const docutils literal notranslate"><span class="pre">'U'</span></code>, only the upper triangular part of <code class="docutils literal notranslate"><span class="pre">A</span></code>
is accessed and the lower triangular part is ignored.</p>
<p>The symbolic factorization is returned as an opaque C object that
can be passed to <a class="reference internal" href="#cvxopt.cholmod.numeric" title="cvxopt.cholmod.numeric"><code class="xref py py-func docutils literal notranslate"><span class="pre">numeric</span></code></a>.</p>
<p>See the comment on
<code class="xref py py-attr docutils literal notranslate"><span class="pre">options['nmethods']</span></code> for details on which ordering is used
by CHOLMOD.</p>
</dd></dl>

<dl class="py function">
<dt id="cvxopt.cholmod.numeric">
<code class="sig-prename descclassname">cvxopt.cholmod.</code><code class="sig-name descname">numeric</code><span class="sig-paren">(</span><em class="sig-param"><span class="n">A</span></em>, <em class="sig-param"><span class="n">F</span></em><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.cholmod.numeric" title="Permalink to this definition">¶</a></dt>
<dd><p>Performs a numeric factorization of a sparse symmetric matrix
as <a class="reference internal" href="#equation-e-chol-ll">(3)</a> or <a class="reference internal" href="#equation-e-chol-ldl">(4)</a>.  The argument <code class="docutils literal notranslate"><span class="pre">F</span></code> is the
symbolic factorization computed by
<a class="reference internal" href="#cvxopt.cholmod.symbolic" title="cvxopt.cholmod.symbolic"><code class="xref py py-func docutils literal notranslate"><span class="pre">symbolic</span></code></a> applied to
the matrix <code class="docutils literal notranslate"><span class="pre">A</span></code>, or to another sparse  matrix with the same sparsity
pattern and typecode, or by
<a class="reference internal" href="#cvxopt.cholmod.numeric" title="cvxopt.cholmod.numeric"><code class="xref py py-func docutils literal notranslate"><span class="pre">numeric</span></code></a> applied to a matrix
with the same sparsity pattern and typecode as <code class="docutils literal notranslate"><span class="pre">A</span></code>.</p>
<p>If <code class="docutils literal notranslate"><span class="pre">F</span></code> was created by a
<a class="reference internal" href="#cvxopt.cholmod.symbolic" title="cvxopt.cholmod.symbolic"><code class="xref py py-func docutils literal notranslate"><span class="pre">symbolic</span></code></a> with <code class="docutils literal notranslate"><span class="pre">uplo</span></code>
equal
to <code class="xref py py-const docutils literal notranslate"><span class="pre">'L'</span></code>, then only the lower triangular part of <code class="docutils literal notranslate"><span class="pre">A</span></code> is
accessed and the upper triangular part is ignored.  If it was created
with <code class="docutils literal notranslate"><span class="pre">uplo</span></code> equal to <code class="xref py py-const docutils literal notranslate"><span class="pre">'U'</span></code>, then only the upper triangular
part of <code class="docutils literal notranslate"><span class="pre">A</span></code> is accessed and the lower triangular part is ignored.</p>
<p>On successful exit, the factorization is stored in <code class="docutils literal notranslate"><span class="pre">F</span></code>.
Raises an <code class="xref py py-exc docutils literal notranslate"><span class="pre">ArithmeticError</span></code> if the factorization does not exist.</p>
</dd></dl>

<dl class="py function">
<dt id="cvxopt.cholmod.solve">
<code class="sig-prename descclassname">cvxopt.cholmod.</code><code class="sig-name descname">solve</code><span class="sig-paren">(</span><em class="sig-param">F</em>, <em class="sig-param">B</em><span class="optional">[</span>, <em class="sig-param">sys = 0</em><span class="optional">]</span><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.cholmod.solve" title="Permalink to this definition">¶</a></dt>
<dd><p>Solves one of the following linear equations where <code class="docutils literal notranslate"><span class="pre">B</span></code> is a dense
matrix and <code class="docutils literal notranslate"><span class="pre">F</span></code> is the numeric factorization <a class="reference internal" href="#equation-e-chol-ll">(3)</a>
or <a class="reference internal" href="#equation-e-chol-ldl">(4)</a> computed by
<a class="reference internal" href="#cvxopt.cholmod.numeric" title="cvxopt.cholmod.numeric"><code class="xref py py-func docutils literal notranslate"><span class="pre">numeric</span></code></a>.
<code class="docutils literal notranslate"><span class="pre">sys</span></code> is an integer with values between 0 and 8.</p>
<table class="docutils align-default">
<colgroup>
<col style="width: 31%" />
<col style="width: 69%" />
</colgroup>
<tbody>
<tr class="row-odd"><td><p><code class="docutils literal notranslate"><span class="pre">sys</span></code></p></td>
<td><p>equation</p></td>
</tr>
<tr class="row-even"><td><p>0</p></td>
<td><p><img class="math" src="_images/math/936adf4ee141541b19f32162898758c0de10fa8d.png" alt="AX = B"/></p></td>
</tr>
<tr class="row-odd"><td><p>1</p></td>
<td><p><img class="math" src="_images/math/36dfce98812a98045ee08f8608788e512eed0fe3.png" alt="LDL^TX = B"/></p></td>
</tr>
<tr class="row-even"><td><p>2</p></td>
<td><p><img class="math" src="_images/math/eb21ee2e4d5e09a543f19ff14ce49d682ad586af.png" alt="LDX = B"/></p></td>
</tr>
<tr class="row-odd"><td><p>3</p></td>
<td><p><img class="math" src="_images/math/7911e1577218dbb02d5e98a311635e1b2347ff57.png" alt="DL^TX=B"/></p></td>
</tr>
<tr class="row-even"><td><p>4</p></td>
<td><p><img class="math" src="_images/math/8621cd51e6bdca11da768ea315ebf1b8c319a975.png" alt="LX=B"/></p></td>
</tr>
<tr class="row-odd"><td><p>5</p></td>
<td><p><img class="math" src="_images/math/e71fcd98134a69534acee44a765e20abefe6e316.png" alt="L^TX=B"/></p></td>
</tr>
<tr class="row-even"><td><p>6</p></td>
<td><p><img class="math" src="_images/math/69ff02759612b6bf7f038a67b9fb4b557f2be247.png" alt="DX=B"/></p></td>
</tr>
<tr class="row-odd"><td><p>7</p></td>
<td><p><img class="math" src="_images/math/38e9e9d3b438f8da362e850c8eb1c0252af8373a.png" alt="P^TX=B"/></p></td>
</tr>
<tr class="row-even"><td><p>8</p></td>
<td><p><img class="math" src="_images/math/00b60df5731428a5a02b49ca7bc4d5aaa7afac9c.png" alt="PX=B"/></p></td>
</tr>
</tbody>
</table>
<p>(If <code class="docutils literal notranslate"><span class="pre">F</span></code> is a Cholesky factorization of the form <a class="reference internal" href="#equation-e-chol-ll">(3)</a>,
<img class="math" src="_images/math/0fcab9067b50b87e868c4fd70f213a086addb964.png" alt="D"/> is an identity matrix in this table.  If <code class="docutils literal notranslate"><span class="pre">A</span></code> is complex,
<img class="math" src="_images/math/249b38285050bda51fd1579d716b5553f8114a0d.png" alt="L^T"/> should be replaced by <img class="math" src="_images/math/183e9814ec579a80697d078963ce0412589322f1.png" alt="L^H"/>.)</p>
<p>The matrix <code class="docutils literal notranslate"><span class="pre">B</span></code> is a dense <code class="xref py py-const docutils literal notranslate"><span class="pre">'d'</span></code> or <code class="xref py py-const docutils literal notranslate"><span class="pre">'z'</span></code> matrix, with
the same type as <code class="docutils literal notranslate"><span class="pre">A</span></code>.  On exit it is overwritten by the solution.</p>
</dd></dl>

<dl class="py function">
<dt id="cvxopt.cholmod.spsolve">
<code class="sig-prename descclassname">cvxopt.cholmod.</code><code class="sig-name descname">spsolve</code><span class="sig-paren">(</span><em class="sig-param">F</em>, <em class="sig-param">B</em><span class="optional">[</span>, <em class="sig-param">sys = 0</em><span class="optional">]</span><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.cholmod.spsolve" title="Permalink to this definition">¶</a></dt>
<dd><p>Similar to <a class="reference internal" href="#cvxopt.cholmod.solve" title="cvxopt.cholmod.solve"><code class="xref py py-func docutils literal notranslate"><span class="pre">solve</span></code></a>, except that <code class="docutils literal notranslate"><span class="pre">B</span></code> is
a class:<cite>spmatrix</cite>, and the solution is returned as an output argument
(as an <code class="xref py py-class docutils literal notranslate"><span class="pre">spmatrix</span></code>).  <code class="docutils literal notranslate"><span class="pre">B</span></code> must have the same typecode as <code class="docutils literal notranslate"><span class="pre">A</span></code>.</p>
</dd></dl>

<p>For the same example as above:</p>
<div class="doctest highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="n">X</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="mi">8</span><span class="p">),</span> <span class="p">(</span><span class="mi">4</span><span class="p">,</span><span class="mi">2</span><span class="p">),</span> <span class="s1">&#39;d&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">F</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">symbolic</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">cholmod</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">F</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">cholmod</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">F</span><span class="p">,</span> <span class="n">X</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="nb">print</span><span class="p">(</span><span class="n">X</span><span class="p">)</span>
<span class="go">[-1.46e-01  4.88e-02]</span>
<span class="go">[ 1.33e+00  4.00e+00]</span>
<span class="go">[ 4.88e-01  1.17e+00]</span>
<span class="go">[ 2.83e+00  7.50e+00]</span>
</pre></div>
</div>
<dl class="py function">
<dt id="cvxopt.cholmod.diag">
<code class="sig-prename descclassname">cvxopt.cholmod.</code><code class="sig-name descname">diag</code><span class="sig-paren">(</span><em class="sig-param"><span class="n">F</span></em><span class="sig-paren">)</span><a class="headerlink" href="#cvxopt.cholmod.diag" title="Permalink to this definition">¶</a></dt>
<dd><p>Returns the diagonal elements of the Cholesky factor <img class="math" src="_images/math/19eef1966f7c545af3ac8c0fa486974d873e3c65.png" alt="L"/>
in <a class="reference internal" href="#equation-e-chol-ll">(3)</a>, as a dense matrix of the same type as <code class="docutils literal notranslate"><span class="pre">A</span></code>.
Note that this only applies to Cholesky factorizations.  The matrix
<img class="math" src="_images/math/0fcab9067b50b87e868c4fd70f213a086addb964.png" alt="D"/> in an <span class="raw-html">LDL<sup><small>T</small></sup></span>
factorization can be retrieved via <a class="reference internal" href="#cvxopt.cholmod.solve" title="cvxopt.cholmod.solve"><code class="xref py py-func docutils literal notranslate"><span class="pre">solve</span></code></a>
with <code class="docutils literal notranslate"><span class="pre">sys</span></code> equal to 6.</p>
</dd></dl>

<p>In the functions listed above, the default values of the control
parameters described in the CHOLMOD user guide are used, except for
<code class="xref c c-data docutils literal notranslate"><span class="pre">Common-&gt;print</span></code> which is set to 0 instead of 3 and
<code class="xref c c-data docutils literal notranslate"><span class="pre">Common-&gt;supernodal</span></code> which is set to 2 instead of 1.
These parameters (and a few others) can be modified by making an
entry in the dictionary <code class="xref py py-attr docutils literal notranslate"><span class="pre">cholmod.options</span></code>.
The meaning of the options <code class="xref py py-attr docutils literal notranslate"><span class="pre">options['supernodal']</span></code>  and
<code class="xref py py-attr docutils literal notranslate"><span class="pre">options['nmethods']</span></code> is summarized as follows (and described
in detail in the CHOLMOD user guide).</p>
<dl class="simple">
<dt><code class="xref py py-attr docutils literal notranslate"><span class="pre">options['supernodal']</span></code></dt><dd><p>If equal to 0, a factorization <a class="reference internal" href="#equation-e-chol-ldl">(4)</a> is computed using a
simplicial algorithm.  If equal to 2, a factorization <a class="reference internal" href="#equation-e-chol-ll">(3)</a>
is computed using a supernodal algorithm.  If equal to 1, the most
efficient of the two factorizations is selected, based on the sparsity
pattern.  Default: 2.</p>
</dd>
<dt><code class="xref py py-attr docutils literal notranslate"><span class="pre">options['nmethods']</span></code></dt><dd><p>The default ordering used by the CHOLMOD is the ordering in the  AMD
library, but depending on the value of <code class="xref py py-attr docutils literal notranslate"><span class="pre">options['nmethods']</span></code>.
other orderings are also considered.
If <code class="docutils literal notranslate"><span class="pre">nmethods</span></code> is equal to 2, the ordering specified
by the user and the AMD ordering are compared, and the best of the two
orderings is used.  If the user does not specify an ordering, the AMD
ordering is used.
If equal to 1, the user must specify an ordering, and the ordering
provided by the user is used.
If equal to 0, all available orderings are compared and the best
ordering is used.  The available orderings include the AMD ordering,
the ordering specified by the user (if any), and possibly other
orderings if they are installed during the CHOLMOD installation.
Default: 0.</p>
</dd>
</dl>
<p>As an example that illustrates <a class="reference internal" href="#cvxopt.cholmod.diag" title="cvxopt.cholmod.diag"><code class="xref py py-func docutils literal notranslate"><span class="pre">diag</span></code></a> and the
use of <code class="xref py py-attr docutils literal notranslate"><span class="pre">cholmod.options</span></code>, we compute the logarithm of the determinant
of the coefficient matrix in <a class="reference internal" href="#equation-e-a-pd">(2)</a> by two methods.</p>
<div class="doctest highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="kn">import</span> <span class="nn">math</span>
<span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">cvxopt.cholmod</span> <span class="kn">import</span> <span class="n">options</span>
<span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">log</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">F</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">symbolic</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">cholmod</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">F</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="nb">print</span><span class="p">(</span><span class="mf">2.0</span> <span class="o">*</span> <span class="nb">sum</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">cholmod</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">F</span><span class="p">))))</span>
<span class="go">5.50533153593</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">options</span><span class="p">[</span><span class="s1">&#39;supernodal&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="mi">0</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">F</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">symbolic</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">cholmod</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">F</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">Di</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">1.0</span><span class="p">,</span> <span class="p">(</span><span class="mi">4</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">cholmod</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">F</span><span class="p">,</span> <span class="n">Di</span><span class="p">,</span> <span class="n">sys</span><span class="o">=</span><span class="mi">6</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="nb">print</span><span class="p">(</span><span class="o">-</span><span class="nb">sum</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">Di</span><span class="p">)))</span>
<span class="go">5.50533153593</span>
</pre></div>
</div>
</div>
<div class="section" id="example-covariance-selection">
<h2>Example: Covariance Selection<a class="headerlink" href="#example-covariance-selection" title="Permalink to this headline">¶</a></h2>
<p>This example illustrates the use of the routines for sparse Cholesky
factorization.  We consider the problem</p>
<div class="math" id="equation-e-covsel">
<p><span class="eqno">(5)<a class="headerlink" href="#equation-e-covsel" title="Permalink to this equation">¶</a></span><img src="_images/math/4e317307fb34d4308336c7711497db123b7d9def.png" alt="\newcommand{\Tr}{\mathop{\bf tr}}
\begin{array}{ll}
    \mbox{minimize} &amp; -\log\det K + \Tr(KY) \\
    \mbox{subject to} &amp; K_{ij}=0,\quad (i,j) \not \in S.
\end{array}"/></p>
</div><p>The optimization variable is a symmetric matrix <img class="math" src="_images/math/52ddc0cde6d632f631533173562fe3ca375b1f32.png" alt="K"/> of order
<img class="math" src="_images/math/5a939c5280da7202ca4531f175a7780ad5e1f80a.png" alt="n"/> and the domain of the problem is the set of positive definite
matrices.  The matrix <img class="math" src="_images/math/7daf0d4815e763eb90f0d5f1dc406f668c1e21db.png" alt="Y"/> and the index set <img class="math" src="_images/math/b988975be41fd13b4d091c10202ba19374643586.png" alt="S"/> are given.
We assume that all the diagonal positions are included in <img class="math" src="_images/math/b988975be41fd13b4d091c10202ba19374643586.png" alt="S"/>.
This problem arises in maximum likelihood estimation of the covariance
matrix of a zero-mean normal distribution, with constraints
that specify that pairs of variables are conditionally independent.</p>
<p>We can express <img class="math" src="_images/math/52ddc0cde6d632f631533173562fe3ca375b1f32.png" alt="K"/> as</p>
<div class="math">
<p><img src="_images/math/e35d7cbe33a34f4456eea34b9bd28642bfc52fa4.png" alt="\newcommand{\diag}{\mathop{\bf diag}}
K(x) = E_1\diag(x)E_2^T+E_2\diag(x)E_1^T"/></p>
</div><p>where <img class="math" src="_images/math/888f7c323ac0341871e867220ae2d76467d74d6e.png" alt="x"/> are the nonzero elements in the lower triangular part of
<img class="math" src="_images/math/52ddc0cde6d632f631533173562fe3ca375b1f32.png" alt="K"/>, with the diagonal elements scaled by 1/2, and</p>
<div class="math">
<p><img src="_images/math/716fbcea39bcb0c1e19a3e26c60836bf5b111b91.png" alt="E_1 = \left[ \begin{array}{cccc}
    e_{i_1} &amp; e_{i_2} &amp; \cdots &amp; e_{i_q} \end{array}\right], \qquad
E_2 = \left[ \begin{array}{cccc}
    e_{j_1} &amp; e_{j_2} &amp; \cdots &amp; e_{j_q} \end{array}\right],"/></p>
</div><p>where (<img class="math" src="_images/math/fc18e7de656eaf3eacbb104912233162cf9042e4.png" alt="i_k"/>, <img class="math" src="_images/math/db5b465a483d35549d1b8a4577cc4163f9ece3aa.png" alt="j_k"/>) are the positions of the nonzero entries
in the lower-triangular part of <img class="math" src="_images/math/52ddc0cde6d632f631533173562fe3ca375b1f32.png" alt="K"/>.  With this notation, we can
solve problem <a class="reference internal" href="#equation-e-covsel">(5)</a> by solving the unconstrained problem</p>
<div class="math">
<p><img src="_images/math/86a2d178f8f888b9c84ff42573fdaf015c8bfbb2.png" alt="\newcommand{\Tr}{\mathop{\bf tr}}
\begin{array}{ll}
\mbox{minimize} &amp; f(x) = -\log\det K(x) + \Tr(K(x)Y).
\end{array}"/></p>
</div><p>The code below implements Newton’s method with a backtracking line search.
The gradient and Hessian of the objective function are given by</p>
<div class="math">
<p><img src="_images/math/a121b0efd1dfa559be59e332b528608d97906c46.png" alt="\newcommand{\diag}{\mathop{\bf diag}}
\begin{split}
\nabla f(x)
    &amp; = 2 \diag( E_1^T (Y - K(x)^{-1}) E_2)) \\
    &amp; = 2\diag(Y_{IJ} - \left(K(x)^{-1}\right)_{IJ}) \\
\nabla^2 f(x)
    &amp; = 2 (E_1^T K(x)^{-1} E_1) \circ (E_2^T K(x)^{-1} E_2)
        + 2 (E_1^T K(x)^{-1} E_2) \circ (E_2^T K(x)^{-1} E_1) \\
    &amp; = 2 \left(K(x)^{-1}\right)_{II} \circ \left(K(x)^{-1}\right)_{JJ}
        + 2 \left(K(x)^{-1}\right)_{IJ} \circ
        \left(K(x)^{-1}\right)_{JI},
\end{split}"/></p>
</div><p>where <img class="math" src="_images/math/ecdff306b10b7e2991a2c20df2fccba0f947155b.png" alt="\circ"/> denotes Hadamard product.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">matrix</span><span class="p">,</span> <span class="n">spmatrix</span><span class="p">,</span> <span class="n">log</span><span class="p">,</span> <span class="n">mul</span><span class="p">,</span> <span class="n">blas</span><span class="p">,</span> <span class="n">lapack</span><span class="p">,</span> <span class="n">amd</span><span class="p">,</span> <span class="n">cholmod</span>

<span class="k">def</span> <span class="nf">covsel</span><span class="p">(</span><span class="n">Y</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Returns the solution of</span>

<span class="sd">         minimize    -log det K + Tr(KY)</span>
<span class="sd">         subject to  K_{ij}=0,  (i,j) not in indices listed in I,J.</span>

<span class="sd">    Y is a symmetric sparse matrix with nonzero diagonal elements.</span>
<span class="sd">    I = Y.I,  J = Y.J.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">I</span><span class="p">,</span> <span class="n">J</span> <span class="o">=</span> <span class="n">Y</span><span class="o">.</span><span class="n">I</span><span class="p">,</span> <span class="n">Y</span><span class="o">.</span><span class="n">J</span>
    <span class="n">n</span><span class="p">,</span> <span class="n">m</span> <span class="o">=</span> <span class="n">Y</span><span class="o">.</span><span class="n">size</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="nb">len</span><span class="p">(</span><span class="n">I</span><span class="p">)</span>
    <span class="n">N</span> <span class="o">=</span> <span class="n">I</span> <span class="o">+</span> <span class="n">J</span><span class="o">*</span><span class="n">n</span>         <span class="c1"># non-zero positions for one-argument indexing</span>
    <span class="n">D</span> <span class="o">=</span> <span class="p">[</span><span class="n">k</span> <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">m</span><span class="p">)</span> <span class="k">if</span> <span class="n">I</span><span class="p">[</span><span class="n">k</span><span class="p">]</span><span class="o">==</span><span class="n">J</span><span class="p">[</span><span class="n">k</span><span class="p">]]</span>  <span class="c1"># position of diagonal elements</span>

    <span class="c1"># starting point: symmetric identity with nonzero pattern I,J</span>
    <span class="n">K</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="n">I</span><span class="p">,</span> <span class="n">J</span><span class="p">)</span>
    <span class="n">K</span><span class="p">[::</span><span class="n">n</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span>

    <span class="c1"># Kn is used in the line search</span>
    <span class="n">Kn</span> <span class="o">=</span> <span class="n">spmatrix</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="n">I</span><span class="p">,</span> <span class="n">J</span><span class="p">)</span>

    <span class="c1"># symbolic factorization of K</span>
    <span class="n">F</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">symbolic</span><span class="p">(</span><span class="n">K</span><span class="p">)</span>

    <span class="c1"># Kinv will be the inverse of K</span>
    <span class="n">Kinv</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="p">(</span><span class="n">n</span><span class="p">,</span><span class="n">n</span><span class="p">))</span>

    <span class="k">for</span> <span class="n">iters</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">100</span><span class="p">):</span>

        <span class="c1"># numeric factorization of K</span>
        <span class="n">cholmod</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">K</span><span class="p">,</span> <span class="n">F</span><span class="p">)</span>
        <span class="n">d</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">F</span><span class="p">)</span>

        <span class="c1"># compute Kinv by solving K*X = I</span>
        <span class="n">Kinv</span><span class="p">[:]</span> <span class="o">=</span> <span class="mf">0.0</span>
        <span class="n">Kinv</span><span class="p">[::</span><span class="n">n</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span>
        <span class="n">cholmod</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">F</span><span class="p">,</span> <span class="n">Kinv</span><span class="p">)</span>

        <span class="c1"># solve Newton system</span>
        <span class="n">grad</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="p">(</span><span class="n">Y</span><span class="o">.</span><span class="n">V</span> <span class="o">-</span> <span class="n">Kinv</span><span class="p">[</span><span class="n">N</span><span class="p">])</span>
        <span class="n">hess</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="p">(</span><span class="n">mul</span><span class="p">(</span><span class="n">Kinv</span><span class="p">[</span><span class="n">I</span><span class="p">,</span><span class="n">J</span><span class="p">],</span><span class="n">Kinv</span><span class="p">[</span><span class="n">J</span><span class="p">,</span><span class="n">I</span><span class="p">])</span> <span class="o">+</span> <span class="n">mul</span><span class="p">(</span><span class="n">Kinv</span><span class="p">[</span><span class="n">I</span><span class="p">,</span><span class="n">I</span><span class="p">],</span><span class="n">Kinv</span><span class="p">[</span><span class="n">J</span><span class="p">,</span><span class="n">J</span><span class="p">]))</span>
        <span class="n">v</span> <span class="o">=</span> <span class="o">-</span><span class="n">grad</span>
        <span class="n">lapack</span><span class="o">.</span><span class="n">posv</span><span class="p">(</span><span class="n">hess</span><span class="p">,</span><span class="n">v</span><span class="p">)</span>

        <span class="c1"># stopping criterion</span>
        <span class="n">sqntdecr</span> <span class="o">=</span> <span class="o">-</span><span class="n">blas</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">grad</span><span class="p">,</span><span class="n">v</span><span class="p">)</span>
        <span class="nb">print</span><span class="p">(</span><span class="s2">&quot;Newton decrement squared:</span><span class="si">%- 7.5e</span><span class="s2">&quot;</span> <span class="o">%</span><span class="n">sqntdecr</span><span class="p">)</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">sqntdecr</span> <span class="o">&lt;</span> <span class="mf">1e-12</span><span class="p">):</span>
            <span class="nb">print</span><span class="p">(</span><span class="s2">&quot;number of iterations: &quot;</span><span class="p">,</span> <span class="n">iters</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
            <span class="k">break</span>

        <span class="c1"># line search</span>
        <span class="n">dx</span> <span class="o">=</span> <span class="o">+</span><span class="n">v</span>
        <span class="n">dx</span><span class="p">[</span><span class="n">D</span><span class="p">]</span> <span class="o">*=</span> <span class="mi">2</span>      <span class="c1"># scale the diagonal elems</span>
        <span class="n">f</span> <span class="o">=</span> <span class="o">-</span><span class="mf">2.0</span> <span class="o">*</span> <span class="nb">sum</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">d</span><span class="p">))</span>    <span class="c1"># f = -log det K</span>
        <span class="n">s</span> <span class="o">=</span> <span class="mi">1</span>
        <span class="k">for</span> <span class="n">lsiter</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">50</span><span class="p">):</span>
            <span class="n">Kn</span><span class="o">.</span><span class="n">V</span> <span class="o">=</span> <span class="n">K</span><span class="o">.</span><span class="n">V</span> <span class="o">+</span> <span class="n">s</span><span class="o">*</span><span class="n">dx</span>
            <span class="k">try</span><span class="p">:</span>
                <span class="n">cholmod</span><span class="o">.</span><span class="n">numeric</span><span class="p">(</span><span class="n">Kn</span><span class="p">,</span> <span class="n">F</span><span class="p">)</span>
            <span class="k">except</span> <span class="ne">ArithmeticError</span><span class="p">:</span>
                <span class="n">s</span> <span class="o">*=</span> <span class="mf">0.5</span>
            <span class="k">else</span><span class="p">:</span>
                <span class="n">d</span> <span class="o">=</span> <span class="n">cholmod</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">F</span><span class="p">)</span>
                <span class="n">fn</span> <span class="o">=</span> <span class="o">-</span><span class="mf">2.0</span> <span class="o">*</span> <span class="nb">sum</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">d</span><span class="p">))</span> <span class="o">+</span> <span class="mi">2</span><span class="o">*</span><span class="n">s</span><span class="o">*</span><span class="n">blas</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">v</span><span class="p">,</span><span class="n">Y</span><span class="o">.</span><span class="n">V</span><span class="p">)</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">fn</span> <span class="o">&lt;</span> <span class="n">f</span> <span class="o">-</span> <span class="mf">0.01</span><span class="o">*</span><span class="n">s</span><span class="o">*</span><span class="n">sqntdecr</span><span class="p">):</span>
                     <span class="k">break</span>
                <span class="n">s</span> <span class="o">*=</span> <span class="mf">0.5</span>

        <span class="n">K</span><span class="o">.</span><span class="n">V</span> <span class="o">=</span> <span class="n">Kn</span><span class="o">.</span><span class="n">V</span>

    <span class="k">return</span> <span class="n">K</span>
</pre></div>
</div>
</div>
</div>


           </div>
           
          </div>
          <footer>
  
    <div class="rst-footer-buttons" role="navigation" aria-label="footer navigation">
      
        <a href="coneprog.html" class="btn btn-neutral float-right" title="Cone Programming" accesskey="n" rel="next">Next <span class="fa fa-arrow-circle-right"></span></a>
      
      
        <a href="fftw.html" class="btn btn-neutral float-left" title="Discrete Transforms" accesskey="p" rel="prev"><span class="fa fa-arrow-circle-left"></span> Previous</a>
      
    </div>
  

  <hr/>

  <div role="contentinfo">
    <p>
        &copy; <a href="copyright.html">Copyright</a> 2004-2020, M.S. Andersen, J. Dahl, L. Vandenberghe
      <span class="lastupdated">
        Last updated on Apr 16, 2020.
      </span>

    </p>
  </div>
  Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/rtfd/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>. 

</footer>

        </div>
      </div>

    </section>

  </div>
  


  <script type="text/javascript">
      jQuery(function () {
          SphinxRtdTheme.Navigation.enable(true);
      });
  </script>

  
  
    
  
    <div class="footer">
      Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 3.0.1.
    </div>

</body>
</html>